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Abstract. An accurate method to compute enclosures of Abelian integrals 
is developed. This allows for an accurate description of the phase portraits 
of planar polynomial systems that are perturbations of Hamiltonian systems. 
As an example, it is applied to the study of bifurcations of limit cycles arising 
from a cubic perturbation of an elliptic Hamiltonian of degree four. 



Nonlinear ordinary differential equations are one of the most common models 
used in any application of mathematical modeUing. In this paper we study famihes 
of such equations 



depending on a small parameter e. 

A fundamental question about such systems is to determine the number and 
location of limit cycles bifurcating from it as e ^ 0. 

In general, the question about the maximal number of limit cycles, and their 
location, of a polynomial planar vector field is the second part of Hilbert's 16th 
problem, which is unsolved even for polynomials of degree 2. For an overview of 
the progress that has been made to solve this problem we refer to [20]. Results for 
the degree 2 case, and a general introduction to the bifurcation theory of planar 
polynomial vector fields can be found in [31]. What is known, is that any given 
polynomial vector field can have only a finite number of Hmit cycles; this is proved 
in [milS]. 

A restricted version of Hilbert's 16th problem, known as the weak, or sometimes 
the tangential, or the infinitesimal, Hilbert's 16th problem, asks for the number of 
limit cycles that can bifurcate from a perturbation of a Hamiltonian system, see 
e.g. [3]. The weak Hilbert's 16th problem has been solved for the degree 2 case, 
see [2]. 

Special cases of Hamiltonian systems are those coming from a one dimensional 
system, H{x, x) — ^ + h{x), which we study in the example given in Section [4Tn If 
one, in addition, assumes that / = 0, and g{x, x) — g{x)x, HI) is known as a Lienard 
equation. Such equations have been thoroughly studied, and the case where dH, 
and g have degree 3 has been solved, see [U [71 [H [9] . We study general g of degree 
3; the set-up of the problem is given in Section 1411 

In this paper we present a rigorous, computer-aided approach to find Hmit cycles 
of planar polynomial vector fields. A different computer-aided approach was intro- 
duced by Malo in his PhD-thesis [22], (also described in [HI US]) which is based on 



1. Introduction 



(1) 



{ 



X = -Hy + ef{x,y) 
y = H^ + eg{x,y), 
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the concept of a rotated vector field, as introduced in [5]. Our approach is com- 
pletely different: we develop a method to rigorously compute what is known as an 
Abelian integral. A brief introduction to Abelian integrals is included in Section [2l 
The concept of a computer-aided proof in analysis is based on techniques to rigor- 
ously enclose the result of a numerical computation. A basis for such a procedure is 
interval analysis, introduced by Moore in [23]. By calculating with sets rather than 
fioating points, it is possible to obtain guaranteed results on a computer, enabhng 
automated proofs for continuous problems. 

We emphasize that the methods developed in this paper are neither restricted 
to any specific degree of the polynomial functions /, and g, nor to the structure 
of the polynomial Hamiltonian H. It can be used to compute AbeHan integrals 
of any polynomial perturbation from any family of compact level curves, ovals, 
of a polynomial Hamiltonian. The method can be used as a computational tool 
to accurately describe the phase portraits of a family of planar systems. In the 
example given in this paper, however, we restrict to the case when / = 0, and dH 
and g have degree 3. 

2. Abelian integrals 

A classical method to prove the existence of limit cycles bifurcating from a 
family of ovals of a Hamiltonian, F/j C H~^{h), depending continuously on h, is to 
study Abelian integrals, or, more generally, the Melnikov function, see e.g. [51 [14]. 
Some caution, however, must be taken regarding the correspondence between limit 
cycles and Abelian integrals, see e.g. [10]. Given a Hamiltonian system and a 
perturbation, 

(2) ( X = -Hy + ff{x,y) 

\ y = + eg{x,y), 

the Abelian integral is defined as 

(3) Hh)^ / f{x,y)dy - g{x,y)dx. 

In this paper all perturbations are polynomial. The most important property of 
Abelian integrals is described by the Poincare-Pontryagin theorem. Let P be the 
return map defined on some section transversal to the ovals of H, parametrised by 
the values h of H, where h is taken from some bounded interval (a, b). We consider 
the displacement function d{h) = P{h) — h. The theorem by Poincare-Pontryagin 
states that 

(4) d{h) = e{I{h) + e(j){h,e)), as e 0, 

where e) is analytic and uniformly bounded on a compact neighbourhood of 

e = 0, h e {a, b). 

3. Computer-aided computation of Abelian integrals 

3.1. Computer-aided proofs. To prove mathematical statements on a computer, 
we need an arithmetic which gives guaranteed results. Many computer-aided proofs, 
including the results in this paper, are based on interval analysis, e.g. |12[ flTl 
[32] . Interval analysis yields rigorous results for continuous problems, taking both 
discretisation and rounding errors into account. For a thorough introduction to 
interval analysis we refer to [T] [23 l [24 l [25 l [27] . 
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3.2. Outline of the approach. The main idea of this paper is to develop a very 
accurate, validated method to enclose the value of a general AbeHan integral. Such 
a method enables us to sample values of I{h). If we can find two ovals T^i, and 
Tfi2 , such that 



then there exists h* G (/ii, /12), such that I{h*) — 0. 

Since P^, the return map of the perturbed vector field, is analytic and non- 
constant, it has isolated fixed points. Thus, we have proved the existence of (at 
least) one limit cycle bifurcating from Th* ■ 

3.3. Computing the integrals. To compute the Abelian integral ([3]) of the form 
(fTBl) . we apply Stokes theorem to get 



where Dh denotes the interior of an oval Th- The reason why we prefer to calculate 
surface integrals, rather than contour integrals, is that we cannot represent the 
ovals of H exactly. We can only find a cover of the ovals, and the area of this cover 
yields the uncertainty of our calculations, automatically handled by the interval 
arithmetic. If we had chosen to compute contour integrals, all of our computations 
would have been subjected to those errors, since we would always integrate over an 
unknown location. When calculating surface integrals, however, the effect of the 
uncertainty of the location of the ovals only contributes on a very small portion of 
the total area of Dh- Note that, inside Dh it is possible to integrate duj exactly, 
that is, there are no truncation errors. 

The actual computation of the integrals is performed in four steps; first we find 
a trapping region for the interesting family of ovals, second we adaptively split this 
region into three parts, one that covers the oval, one representing the inside and 
one representing the outside, third we change the coordinates on the boxes covering 
the oval in order to minimise the area of the cover, fourth we integrate duj on the 
boxes representing the inside and the cover of the oval. 

The first step is simple, since we only consider ovals that are situated inside a 
homo- or heteroclinic orbit. A short branch-and-bound algorithm quickly finds a 
box enclosing the homo- or heterochnic orbit; this box is our initial domain used 
for the main part of the program. 

In the second step - the adaptive splitting of the domain - we perform a series 
of tests to determine whether a box B intersects the oval, is inside it, or outside 
it. We start by evaluating the Hamiltonian on B using monotonicity and central 
forms; since the Hamiltonians we study are sufficiently simple, we implement the 
derivatives symbohcally. Three cases occur: if H < h, then B is inside the oval, and 
we label B as such. \{ H > h then B is outside the oval and we ignore it. Finally, 
if h G H, then we try to perform the change of variables as described below. If the 
change of variables procedure fails, and the size of B is greater than some stopping 
tolerance, minsize, then we split the box B into four parts and re-examine them 
separately. If the size of B is smaller than minsize, then it is labelled fail. If the 
change of variables procedure works, then we label B as on. 



(5) 



I{hi)I{h2) < 0, 



(6) 




dw 
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Figure 1. The labelling of boxes intersecting an oval. 

The third, and most complicated, part of our program is the change of variables 
in the boxes that intersect the oval. Let b & B he the midpoint of B. Compute 

u = yH{h), 

and choose v such that 

u _L w and vi > 0. 

Using the labelling illustrated in Figure [TJ let right and left be the sides in- 
tersected by the line b + tv, t e M.. Denote the intersection points by p, and q, 
respectively. The possible configurations of an intersection of the oval with a box 
are illustrated in Figure [2l The restriction of H to the sides right and left, re- 
spectively, are one-dimensional functions, and the location of the intersections can 
be approximated, and their uniqueness proved, using the interval Newton method 
[23] initiaHsed from the points p, and q, respectively. Let, 

accuracy = minsize/10. 

Define the points Pup, Pdown on the right-side and the points g„p, qdown on the 
left-side at the distance accuracy from p and q, respectively, as illustrated in 
Figure [31 




Figure 2. The possible configurations of the intersection of an 
oval and a box. 



If the following conditions hold, then the oval is inside the tube illustrated in 
Figure O and we can change coordinates to get a small box, which is guaranteed to 
contain the segment of the oval passing through B. This small box represents the 
error caused by the unknown location of the oval. 
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Figure 3. Constructing a small, local enclosure of the oval. 
Condition 3.1. 

sign {H{pup) - h) = sign(i/((j„p) - h) 

= -sigIl{H (pdown) - h) . 

= -sign{H{qdown) - h) 

Let lup'i and Idown denote the line segments between pup and qup, and pdown and 
qdown, respectively. Denote by H' differentiation with respect to the parametrisa- 
tion of the line lup, and Idown, respectively. 

Condition 3.2. 

^ (H{lup) - h) or ^ H'iLp) 
i {H{ldou,n) -h) or ^ H'{1 

down ) • 

Let othersidel, and otherside2 be the two other sides of the box B, that is, 
othersidel U otherside2 U right U left = {1, 2, 3, 4}. 
Condition 3.3. 

Th n othersidel = and Th H otherside2 = 0, 

We enclose the segment of the oval inside of the box between two straight lines: 
Condition 13 . II guarantees that the points Pup, Pdown, lup, and qdown are on different 
sides of the oval as in Figure O Condition 13.21 guarantees that the lines lup and 
Idown do not intersect the oval, and Condition [33] guarantees that the oval does not 
cross the other sides of the box. Recall that the uniqueness of p and q is proved 
as they are approximated. Hence, we have proved that the segment of the oval 
crossing the box has exactly two intersections with the boundary of the box, and 
that it is confined to the region between lup and Idown- 

If (|3.ip . I|3.2p . and (|3.3p . hold, then we set accuracy=accuracy/2, re-calculate 
Pup: Pdown, qup, and qdown, and try to verify l|3.ip . p.2p . and p.3p . This procedure 
is iterated until (|3.ip . or l|3.2p do not hold. Finally, we label B as on. 

The fourth and final part of our integration algorithm, is the actual integration. 
The integration is done separately for the boxes that are labelled, inside, fail, 
and on. 

If B is inside we compute 

'sup(Bi)^+i inf(Bi)'+i\ /sup(B2)^+i inf(B2p+^~ 



(7) / x'w^ dxAdy = , , . , 

If B is labelled fail, we know that B might intersect the oval, that is, we have 
neither been able to prove intersection, nor non-intersection. Therefore, we must 
include any possible result; the integral over B is calculated as the interval hull of 
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and the largest, and smallest, respectively, result of ^ calculated on a subbox 
B C B. 

Boxes labelled fail cause large over-estimations. Fortunately such boxes are 
rare, typically less than 5% of the on-boxes, see Section |4l If minsize is taken 
sufficiently small, the effect of the fail-boxes is negligible. 




Bi 



Figure 4. The change of variables splitting. 



The boxes that are labelled on, are split into five parts, as illustrated in FigurelH 
By construction, none of the triangles, Ti,Tu, or boxes Bi, Bu in the splitting of B 
intersect the oval, thus it suffices to evaluate H in one point of each, and hence they 
can all be labelled as inside or outside. The boxes Bi,Bu are then treated as 
above, that is, if they are labelled inside they are integrated according to ([7]), and 
if they are labelled outside they are neglected. A triangle labelled outside is also 
neglected, the integrals on triangles labelled inside are enclosed by the formula 

(8) / x'y^ dxAdy eDTlDTi\T\, 

Jt 

where DT is the box hull of T, and |T| is the area of T. This gives a reasonably 
narrow enclosure of the integral, since the width of DT is typically small. The 
parallelepiped, P, which covers the segment of the oval, remains to be studied. 
When we integrate over P, the same problem as in the fail case occurs; we do not 
know how much of the parallelepiped to include. Therefore, we have to take the 
hull of all possible outcomes. Hence, the integrals are computed as 

(9) j x'y^ dxAdye Hull (o, nPinPl\P\J , 

where DP is the box hull of P and \P\ is the area of P. 

The value of the Abelian integral is enclosed by summing over all the computed 
integrals that are labelled as either inside, fail, or on. 

ide *^ 

(10) + Esefaii Hull(0, 0) 

4. Computational results 

In this section we apply the methods developed in Section 13.31 to an eUiptical 
Hamiltonian of degree four, described in Section 14.11 The main idea is to inte- 
grate monomial forms at some points, and then to specify the coefficients of the 
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perturbation lo such that I{h) is zero at the sampled points. Therefore, let 
(11) I,j{h) = I x'y^ dxAdy. 



We sample at some number of /i-values, uniformly distributed between the saddle 
loops and the singularity. From these calculations we deduce candidate coefficients. 

Given some candidate coefficients of the form uj, we calculate the Iij{h), at 
intermediate ovals. If the linear combination of the Iij{h) has validated sign changes 
between the sample points we are done: it has been proved that the corresponding 
perturbation yields bifurcations with the given number of limit cycles as e — »■ 0. 

All computations were performed on a Intel Xeon 2.0 Ghz, 64bit processor with 
7970Mb of RAM. The program was compiled with gcc, version 3.4.6. The software 
for interval arithmetic was provided by the C-XSC package, version 2.1.1, see .18|. 

4.1. Example - bifurcations from a figure eight loop. We study the elliptic 
Hamiltonian of degree 4 with a figure eight loop, given by 

. , 11^ x^ 1 — A q A , 

(12) H^^ + - + ^x^--x\ 

where A e (0, 1), see [9]. The corresponding differential system has two centres, at 
H — —j^{2X + 1), and H — —j^X^{X + 2), that are surrounded by a figure eight 
loop, located at H = 0, see Figure [H As A grows the right loop grows; A = 1 is a 
symmetric figure eight loop. We choose to study A = 0.95; a motivation why we 
want A large is given below. The Hamiltonian (fT2| corresponds to the differential 
system, 

X = -Hy = -y 



(13) 



y = = x'-^ + {1- X)x^ - Xx 



We are interested in Hmit cycles bifurcating from the periodic solutions of (fT3|) . 
corresponding to integral curves of lfT2|) . The closed level-curves of (fT2l) are called 
ovals. In a series of papers [UlTlIHlI^, Dumortier and Li study cubic perturbations 
of elliptic Hamiltonians corresponding to Lienard equations. That is, 

(14) X + e{a + f3x + ^x^)x + ax^ + bx^ + ex = 0. 

For the elliptic Hamiltonians of degree four with compact ovals, there are five 
different classes of phase portraits, see e.g [3]. They are, the truncated pendulum, 

the saddle loop, the global centre, the cuspidal loop, and the figure-eight loop. 

3 

Compared to the Lienard case, we add a fourth term, 5^, to the perturbation, 
and explore what kind of bifurcations we can prove to exist. We study the perturbed 
system. 



(15) 



X ^ -y 



y = x^ + {1 - X)x^ - Xx + e (^{a + Px + ^x^)y + S^'j . 
The 1-form associated with this perturbation is 

(16) uj = - (^{a + (3x + 7x^)2/ + (5y^ dx. 

For computational efficiency we primarily study its exterior derivative, 

(17) du) = {{a + l3x + jx'^) + Sy^) dx /\dy. 
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Figure 5. The elliptic Hamiltonian of degree 4 with a figure eight 
loop, studied in section \4A\ 



In [29] Petrov proves that when restricting to one family of ovals, surrounding 
one of the two centres, the space of Abelian integrals has dimension 4, and that 
the space has the Chebyshev property, that is, the number of zeros of a function in 
this space is less than the dimension of the space. He also proves that this bound 
is sharp. To construct an example with more than three limit cycles surrounding 
either of the two centres, we can therefore not simply use the Chebyshev property 
of the space of Abelian integrals. 

Our heuristic argument to guess parameters is the following: we start by in- 
tegrating at 100 uniformly distributed ovals, in each eye of the loop. We do this 
with moderate accuracy, which gives a fast and sufficiently precise result. Since 
we have chosen to study a figure eight loop that is not far from being symmetric, 
it is reasonable to assume that the two branches behave similarly, which makes it 
probable that we should be able to determine coefficients so that each branch has 
two zeros. To determine such zeros, we solve the following linear system: 

1 

-1 
1 

-1 





-0.0362) 


4o( 




-0.1208) 


Hoi 


lU 


-0.1812) 


Hoi 




-0.1054) 


Hoi 



0.0362) 72o(--0.0362) 

0.1208) 7^0 (-0.1208) 

0.1812) 4o(-0.1812) 

0.1054) JJo (-0.1054) 



/o2 (-0.0362) 1 r Q 

/^2(-0.1208) 13 

/^2(-0.1812) 7 

7^2 (-0.1054) J [ S 



where llj{h), and Ilj{h), denote the monomial Abelian integrals calculated on the 
left and right ovals, respectively. 

This gives the approximate solution a = 438.4905, /3 = -25.2469, 7 = -452.7899, 
and b = —741.0341, which we use as our perturbation. The graph of the resulting 
function is given in Figure [7l which appears to have 4 zeros, illustrated in Figure 
[6l This, of course, has to be proved. 

To prove that the perturbation constructed above has 4 zeros, we procede as in 
the previous examples, and compute enclosures of the Abelian integral at intermedi- 
ate ovals. On the left branch we calculate /(-0.0121), /(-0.0846), and /(-0.1933), 
and on the right branch we compute /(-0.0105), /(-0.0738), and /(-0.1686). The 
result is given in Tables [U and [21 
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~-?.5 -1 -0.5 0.5 1 1.5 

Figure 6. The limit cycles from Example 14. Ij unstable limit cy- 
cles are dashed. Note that we only prove the existence of the limit 
cycles, their locations as drawn in the figure are the locations of 
the ovals they bifurcate from. 

20p 
15- 
10- 
5- 



0- 



-!^25 -0.2 -0.15 -0.1 -0.05 

Figure 7. The two branches of the Abelian Integral for the figure 
eight loop. 

at the intermediate ovals, 

[+8.698, +9.290], 
[-2.204,-1.780], 
[+0.9121, +1.119], 
[+11.56, +12.10], ■ 
[-1.181,-0.7959], 
[+0.2095, +0.3847] 

Hence, the system with the given perturbation has four limit cycles, one attracting 
and one repelling inside each loop, see Figures [6] and [8l The run-time of the program 
was, for the left (right) branch, 82 (78) seconds, a total of 1182 (1166) boxes were 
used to cover the 3 ovals, 82 (56) of these belong to the fail class. 

To prove that the unstable separatrices of the saddle are attracted to a limit 
cycle enclosing the figure eight loop, as indicated in Figure [3 we first calculate 
I°{h), the outer Abelian integral, for some h > values with low accuracy to find 





Finally, we compute I (h), and I^{h) 

/'(-0.0121) = 

/' (-0.0846) = 

/'(-0.1933) = 

^ ' /'■(-0.0105) = 

/'■(-0.0738) = 

r(-o.u 
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h 


^00 






^02 


-0.0121 


[1.206,1.207] 


[-1. 034,-1. 033[ 


10.9945,0.9951] 


10.1290,0.1293] 


-0.0846 


[0.7661,0.7665] 


[-0. 7073,-0. 7068[ 


10.6902,0.6907] 


10.05829,0.05839] 


-0.1933 


10.2219,0.2222] 


[-0.2178,-0.2175[ 


10.2160,0.2164] 


10.005326,0.005340] 


Table 1 . The computed enclosures for the left branch of the figure 


eight loop. 








h 


rr 

-'oo 


ho 


rr 
^20 


rr 

-"02 


-0.0105 
-0.0738 
-0.1686 


]1. 077,1. 078] 
[0.6846,0.6851[ 
[0. 1984,0. 1987[ 


[0.8773,0. 8778[ 
[0.6002,0.6006] 
10.1848,0.1850] 


10.8033,0.80.39] 
10.5573,0.5577] 
10.1744,0.1747] 


10.1006,0.1008] 

10.04545,0.04553] 

10.004154,0.004164] 



Table 2. The computed enclosures for the right branch of the 
figure eight loop. 



h 


-'oo 


r° 

ho 


T° 

-'20 


^02 


0.09 
0.11 


13.576,3.587] 
13.776,3.786] 


1-0.1843,-0.1709] 
1-0.1862,-0.1740] 


12.560,2.575] 
12.708,2.724] 


10.5307,0.5376] 
10.6044,0.6109] 



Table 3. The computed enclosures for the outside of the figure 
eight loop. 



an indication of a sign change. It appears that a limit cycle bifurcates from an oval 
close to H = 0.1. Therefore, we compute /°(0.09), and /°(0.11), 

/°(0.09) = [+8.715, +24.83], 
/°(0.11) = [-25.37,-9.821]. 

These calculations verify that the perturbed system has an attracting Hmit cycle 
bifurcating from an oval outside the figure eight loop. The run-time of the program 
was 39 seconds, a total of 496 boxes were used to cover the 2 ovals, 52 of these 
belong to the fail class. 

5. Conclusions 

We have presented a method to rigorously calculate Abelian integrals. The 
method can be applied to study any polynomial perturbation of a polynomial 
Hamiltonian vector field. As an application, we have applied the method to an 
elliptic Hamiltonian of degree 4. 

The method can be used in several ways: either one can use it to verify that a 
specific perturbation guessed by some other method indeed has a certain number of 
zeros, or one can use it as in Section |4!T] to sample and plot the monomial Abelian 
integrals. In the latter case, if a good choice of parameters can be made from the 
approximate knowledge of the monomial AbeHan integrals, then one can re-use the 
program to verify that guess, as is done in Section 

A major challenge is to device a method which can be used to guess what per- 
turbations to investigate. One such method that appears in the literature is that of 
a detection function, as used in e.g. [S^- Another problem, which we have ignored 
in this paper, is that typically when one has a Hamiltonian depending on param- 
eters, the maximal number of limit cycles that can bifurcate from one member of 
this family, will only appear for some special values of the parameters. It would 



COMPUTER-AIDED COMPUTATION OF ABELIAN INTEGRALS 



11 




Figure 8. The perturbed figure eight loop, here with e = 0.001, 
illustrating the 5 limit cycles found in Section 14.11 

therefore be desirable to develop conditions indicating how to choose one candidate 
system from a family. In Section 14.11 we give a completely heuristic argument why 
we want to have A large. 
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